Direct on-swab metabolic profiling of vaginal microbiome host interactions during pregnancy and preterm birth

The pregnancy vaginal microbiome contributes to risk of preterm birth, the primary cause of death in children under 5 years of age. Here we describe direct on-swab metabolic profiling by Desorption Electrospray Ionization Mass Spectrometry (DESI-MS) for sample preparation-free characterisation of the cervicovaginal metabolome in two independent pregnancy cohorts (VMET, n = 160; 455 swabs; VMET II, n = 205; 573 swabs). By integrating metataxonomics and immune profiling data from matched samples, we show that specific metabolome signatures can be used to robustly predict simultaneously both the composition of the vaginal microbiome and host inflammatory status. In these patients, vaginal microbiota instability and innate immune activation, as predicted using DESI-MS, associated with preterm birth, including in women receiving cervical cerclage for preterm birth prevention. These findings highlight direct on-swab metabolic profiling by DESI-MS as an innovative approach for preterm birth risk stratification through rapid assessment of vaginal microbiota-host dynamics.

T he vaginal microbiome is a key mediator of reproductive tract pathophysiology. Unlike other mucosal surfaces such as the gut, low diversity in the vaginal microbiome and dominance by Lactobacillus species is considered a hallmark of health, particularly during reproductive years. Lactobacillus species depletion and increased microbial diversity are characteristic of bacterial vaginosis (BV) 1 , and associates with both the increased risk of acquisition 2 and ineffective treatment of sexually transmitted infections, including HIV 3,4 . High diversity, BVassociated vaginal microbiota have also been linked to HPV infection and cervical dysplasia 5,6 , in vitro fertilization failure 7 , miscarriage 8 and preterm birth [9][10][11][12][13] . Current clinical diagnosis of vaginal infection is largely limited to subjective assessment of clinical symptoms in addition to time-consuming microscopic evaluation of vaginal swab samples, pH and, in some cases, culture [14][15][16] . These approaches fail to detect many clinically relevant species and lack insight into overall community composition. Such information can be obtained using next-generation sequencing-based approaches (e.g. metagenomics and metataxonomics), but these involve complex and time-consuming sample preparation, data processing and expense that prevents their use for routine bedside testing. Further, molecular-based characterization of microbiota is unable to assess microbiota:host interactions that ultimately determine health and disease phenotypes. Therefore, there exists a need for rapid, point-of-care vaginal diagnostics to facilitate faster clinical decision making, more judicious use of antibiotics and targeted treatment strategies.
A common mechanistic pathway linking sub-optimal vaginal microbiota composition (VMC) and pathophysiology is activation of host-innate immune response and inflammation 17 , which can be suppressed by Lactobacillus species such as L. crispatus, through the modification of the metabolic milieu of the cervicovaginal mucosa 16,18 . During pregnancy, untimely activation of inflammation in gestational tissues (e.g. cervix, foetal membranes) caused by ascending vaginal infection is thought to cause a significant proportion of infection-associated preterm birth 19,20 , which remains the primary cause of death in children under 5 years of age 21 . Consistent with this, we have recently reported that cervical cerclage, a procedure used to reinforce the cervix in women at risk of preterm delivery due to cervical shortening, can induce vaginal bacterial dysbiosis, inflammatory activation and premature cervical ripening associated with increased risk of preterm birth, if performed using a braided instead of monofilament cerclage material 22 . In contrast, L. crispatus dominance of the vaginal niche during pregnancy has been associated with protection against preterm birth 10,12,22 . Despite this, individual and ethnic variation in VMC and host response 11,23,24 as well as cost and time constraints have limited the utility of current analytical methods used for vaginal microbiota characterization to inform clinical decision-making during pregnancy. Current methods for preterm birth prediction have poor positive predictive value and fail to provide insight into underlying aetiology, which may explain low efficacy of interventions designed to prevent preterm birth 25 .
We have recently described a method that enables rapid, objective assessment of the chemical composition of mucosal surfaces using sample preparation-free, direct on-swab desorption electrospray ionization mass spectrometry (DESI-MS) 26,27 . This involves directing a pneumatically assisted electrospray of charged aqueous droplets directly onto a rotating swab, where it forms a liquid film that desorbs and ionizes molecules from the sample, which are transferred to a mass spectrometer via an atmospheric pressure ion-transfer line. With this method, metabolic profiles can be acquired rapidly (in <3 min) and directly from a swab sample without the need for laborious sample preparation or chemical extraction procedures. While lack of chromatographic separation limits assay selectivity, it provides the advantage of vastly simplifying the amount of maintenance and operator intervention required. Additionally, since the solvent stream only ablates sample from a small area of the swab, the method is virtually non-destructible supporting multi-assay use from the sampling device (e.g., cytokine/immune marker profiling or 16S rRNA sequencing) 16 . Collectively, these characteristics make direct on-swab DESI-MS particularly suitable for deployment in clinical point-of-care settings.
We hypothesized that vaginal microbiota:host interactions would be reflected in the metabolic milieu of the cervicovaginal mucosa during pregnancy and provide sufficient biochemical information to predict both bacterial composition and host immune response. To test this, direct on-swab DESI-MS profiling was used to characterize the cervicovaginal metabolome of two independent cohorts of women prospectively sampled throughout pregnancy (VMET, n = 160; 455 swabs; VMET II, n = 205; 573 swabs, Fig. 1A). These data were then integrated with matched metataxonomic and immuno-profiling data to identify DESI-MS metabolic signatures predictive of VMC and local inflammatory status. We then examined if these predictive models could be used to monitor changing VMC and host inflammatory responses that associate with preterm birth and clinical interventions (e.g., cervical cerclage) used during pregnancy.

Results
Baseline characteristics of the prospective study subjects. A total of 160 women recruited to the VMET study were longitudinally sampled throughout pregnancy (n = 455 swabs) (Fig. 1A). These samples were all subjected to metabolic profiling by DESI-MS and a cross-platform comparison using five complementary liquid chromatography-MS (LC-MS) assays (Fig. 1A). A second independent patient cohort (n = 205, VMET2 study) were longitudinally sampled at comparable time points (n = 573 swabs). Samples from these patients were analysed by DESI-MS and a subset were used for immuno-profiling studies. Key demographics including maternal age, BMI, gravida and ethnicity were similar between the two patient cohorts (Table 1). Preterm births rates in these high-risk populations were 18% (n = 29; VMET) and 21% (n = 44; VMET2).
Using linear mixed effect modelling, a total of 113 metabolite features determined in DESI-MS negative and positive ion modes were found to significantly discriminate between LDOM and LDEPL in independent analyses of both patient cohorts (Fig. 1B) Fig. 1 Direct on-swab desorption electrospray ionization mass spectrometry (DESI-MS) metabolic profiling of cervicovaginal fluid enables robust prediction of vaginal microbiome compositions. A Study design and longitudinal multi-omic sampling and analysis workflow of cervicovaginal swab samples collected from two independent pregnancy cohorts (VMET: n = 165, 455 swabs; VMET2: n = 205, 573 swabs). Data from each cohort were analysed independently, and features selected only if their Benjamini-Hochberg q-value was smaller than 0.05 in both datasets, after matching the hits from each analysis by their m/z values. B Heatmaps representing relative concentrations of DESI-MS (negative mode)-derived metabolic features (n = 88) significantly differing between Lactobacillus spp. dominated (L-dominated, green) and Lactobacillus spp. depleted (L-depleted, red) states in both independent patient cohorts (see Table S2). C Boxplots of representative discriminatory metabolic features with corresponding Benjamini-Hochberg q-values identified including thiomalic acid, leucyl-serine, docosanoic acid (C22:0), lignoceric acid (C24:0) with calculated z-score measured in the two patient cohorts VMET2 (left, n = 203, 539 swabs) and VMET (right, n = 160, 428 swabs). The lower and upper bounds of the box represent the 25th and 75th percentile values, respectively, and the interior horizontal line the median value. Whiskers are drawn from the corresponding box boundary to the most extreme data point located within the box bound ± 1.5 × IQR (interquartile range). m/z mass-to-charge ratio, CST community state type.  Fig. 3). In addition, the proportion of total metabolome variance explained by each factor was estimated with permutational multivariate analysis of variance (PERMANOVA, Supplementary Table 2). While the R 2 values estimated for each covariate varied between metabolic profiling assays, their relative importance was consistent, with between-individual variability explaining the most variance (35-45.2%), followed by CST (2.2-25.2%), ethnicity (0.5-2.4%) and gestational age (0.4-1.8%). Of those, thiomalic acid and leucyl-serine were consistently higher in LDOM samples, whereas docosanoic acid and lignoceric acid were significantly higher in LDEPL samples (Fig. 1C). Random forest classifier and ROCcurve analysis were then used to assess the performance of DESI-MS to predict VMC ( Fig. 2 and Supplementary Fig. 4). Robust prediction of genera-level classification (LDOM vs LDEPL) was observed across patient cohorts particularly using features derived from negative ion polarity mode (VMET/VMET2; AUC 94. In vitro DESI-MS profiling of vaginal commensal and pathogenic bacterial isolates. We next investigated if any of the discriminatory metabolites identified in the in vivo DESI-MS analyses could also be observed in culture isolates (n = 25) of bacterial species recognized as being predominant members of major vaginal CSTs. In total, 27 of the discriminatory metabolites identified in vivo were detected by DESI-MS in swabs of culture biomasses following correction for background media concentrations (Fig. 3). Approximately half of these metabolites were detected at levels lower than that observed in media background controls, whereas the remainder were found at levels higher than media background. No clear relationship was observed between the DESI-MS detected levels of these metabolites in vitro, CST membership nor levels observed in vivo.
Assessment of host immune response at the mucosal interface using direct on-swab DESI-MS. A panel of 22 soluble immune markers including cytokines, chemokines, immunoglobulins and members of the complement system were measured in a subset of the VMET2 cohort samples (n = 391). Random forest regression analysis was used to predict the log-transformed concentrations of each marker using DESI-MS-derived features. Robust prediction (cross-validated R 2 > 0.25) was observed particularly for IL- 4A). Immune markers with CV R 2 > 0.1 (n = 9) and DESI-MS features with an R 2 > 0.1 for the linear regression against the immune marker (n = 23) revealed a metabolic signature strongly associated with local immune phenotype primarily characterized by altered levels of long-chain fatty acids, glycerophospholipids and ceramides ( Fig. 4B and Supplementary Table 4). Both DESI-MS predicted and immunoassay-measured levels of C3b, IL-1β, IgG2 and IgG3 were elevated in Lactobacillus-depleted vaginal microbiomes indicating activation of the local innate and adaptive immune response (Fig. 4C).
Direct-on swab DESI-MS as a potential tool for PTB risk stratification. We next examined the relationship between pregnancy outcome, VMC and inflammatory status as predicted by DESI-MS by combining all available data from both VMET and VMET2 patient cohorts. High vaginal microbiome diversity and instability during pregnancy, as defined by shifts between Lactobacillus-dominated and Lactobacillus-depleted compositions classified with 16S rRNA gene-based metataxonomics, was associated with an increased risk of preterm birth compared to those women maintaining Lactobacillus-dominance throughout pregnancy (odds ratio (OR) 1.97, 95% confidence interval (CI) 1.03-1.66, p = 0.04) (Fig. 5A). Similarly, vaginal microbiota instability determined solely by DESI-MS was also associated with a higher risk of preterm birth although this did not reach statistical significance (OR 1.47, 95% CI 0.75-2.78, p = 0.25).
Increased levels of IL-1β, measured directly by immunoassay or predicted by DESI-MS, were associated with high-diversity VMC in both those women experiencing term birth or preterm birth. However, levels were significantly higher in preterm birth cases (preterm 59.63 pg/ml vs term 1.94 pg/ml, p = 0.006, Welch two sample t-test, Fig. 5B). Increased levels of MBL, again either measured directly by immunoassay or predicted by DESI-MS, were also observed in high-diversity VMC of women  Table 3 and Supplementary Fig. 4).
We next tested if DESI-MS metabolic, metataxonomics and inflammatory marker profiles obtained at three different stages of gestation, could predict subsequent preterm birth. Overall, the predictive capacity of these models was poor (Supplementary Table 5). We therefore focused our subsequent analyses on preterm birth phenotypes more likely to associate with dysregulated vaginal microbiota-host interactions, as previously described 22 . In women at high risk of preterm birth due to cervical shortening, vaginal levels of MBL were most frequently increased in those treated with a cervical cerclage using braided suture material compared to those where cerclage had been performed with monofilament suture (10/11, 91% vs 9/21, 43%, p = 0.011, Fisher's exact test). This was similarly detected when VMC and MBL levels were estimated using only DESI-MS analysis of vaginal swabs (7/11, 64% vs 8/21, 38%, p = 0.316; Fig. 5D). Immunoassay also indicated increased IL-1β after braided cerclage insertion, but this was less consistently detected by DESI-MS prediction (Fig. 5E). However, following treatment with braided cerclage material, DESI-MS did accurately predict low levels of IL-1β in women delivering at term compared to those who subsequently delivered preterm (Fig. 5F).

Discussion
Despite recent developments in metataxonomics and metagenomics, VMC characterization in clinical settings remains largely limited to culture and microscopy, which like molecular-based approaches, fails to capture information regarding host response. Our method, which is easily amenable to bedside point-of-care testing 31,32 , addresses this limitation by leveraging information contained within the cervicovaginal mucosa to provide robust detection of VMCs and simultaneous estimation of host immune and inflammatory status.
Many of the discriminatory metabolites identified in our study have been previously reported using gas chromatography-MS or LC-MS based assays in women suffering from BV [33][34][35][36] or in HPV infection 37 and have biologically plausible roles in mediating vaginal health and disease. For example, increased vaginal levels of short-and long-chain fatty acids and biogenic amines are associated with activation of pro-inflammatory pathways 36,38-41 , which contribute to reduced barrier integrity of the epithelia and consequently, increased risk of infection 34,[42][43][44][45] . Furthermore, cervicovaginal lipid species have recently been associated with tumour progression and genital inflammation defined using an aggregated score, in a study of 78 non-pregnant, HPV-negative/ positive women with cervical dysplasia 37 . Lipid changes in vaginal discharge have also been associated with vulvovaginal candidiasis 46 . In our study, we profiled a larger panel of soluble immune mediators, including humoral response mediators, and identified associations with the metabolome via regression models. Despite different experimental approaches, our results provide further evidence that lipid species are a core component of the metabolic signature of inflammation and immune response in the vaginal niche. Some of the discriminatory markers used in our in vivo DESI-MS prediction models were also detected by DESI-MS in vitro swab analyses of bacterial isolate biomasses. It should be noted that the bacterial isolates analysed here were not representative of the genetic diversity of vaginal commensal and pathogenic strains; however, our analyses indicate that the metabolite signatures predictive of microbiota composition and inflammation in vivo are likely derived from both bacterial and host sources.
Some studies have estimated that around 23-30% of total cervicovaginal metabolic variation is associated with bacterial composition 36,47,48 . Similar estimates obtained in our data emphasize the impact of the microbial composition on the metabolome, which is greater than the effect of gestational age, ethnicity, maternal age, or BMI. However, using the repeated measures per individual, we found that the top proportion of variance explained could be attributed to between-individual variability. The magnitude of the R 2 values for microbial composition was higher in the LC-MS assays compared to DESI-MS. This may represent greater coverage of relevant constituents of the metabolome that are not detected by DESI-MS or may reflect the technical variability of DESI-MS which, like other ambient ionization methods, is recognized to be larger than that of LC-MS where implementation of experimental quality control (QC), instrumental drift correction and data filtering procedures are more established 49 . Despite this, our findings highlight the capacity of direct on-swab DESI-MS to rapidly capture this information from the metabolome without the need for laborious sample preparation and comparatively high per/assay time and economic costs associated with coupled chromatography-MS assays for similar applications. This lends itself to customization and automation for ease-of-use, which are important characteristics for point-of-care testing 32 . Our approach also offers the advantage of providing objective simultaneous assessment of microbiota composition and inflammatory state. In comparison, current 'gold standard' diagnostic methods for vaginal infection (e.g. Amsel criteria) are limited to a combined subjective assessment of clinical symptoms and microscopy grading of bacterial morphotypes as well as selective culture.
There is now substantial data supporting a role for the vaginal microbiome and host immune responses in shaping preterm birth in a proportion of women [10][11][12]22 . However, prediction of preterm birth using DESI-MS metabolic profiles, metataxonomics or inflammatory marker data in our patient cohort was poor. This is not surprising given the fact that preterm birth is a multiaetiological disease state that can be caused by many different factors, including non-microbial and non-immune related causes 19,20 . Because of this, we focused subsequent analyses on women who receive cervical cerclage with braided suture, who we have previously shown are at increased risk of preterm birth that involves a phenotype characterized by vaginal dysbiosis and local immune activation 22 . Here we confirmed these findings in an independent patient population and highlight the utility of direct swab analysis by DESI-MS to rapidly detect both microbiota and inflammatory changes caused by the intervention that associate with subsequent preterm birth risk. The ability to provide such information at point-of-care would be transformative for directing clinical decision making and ultimately improving outcomes for these women and their babies.
In this study, DESI-MS prediction of VMCs was capable of monitoring vaginal microbiome diversity and instability longitudinally throughout pregnancy, which in the study participants was associated with increased risk of preterm birth. These results are consistent with a recent meta-analysis of metataxonomicsbased studies of the vaginal microbiome in pregnancy, which reported higher variance of VMC across trimesters in women subsequently experiencing preterm birth 50 . An additional strength of our approach is the ability to simultaneously capture relationships between VMC and local immune response. Innate immune activation in the vagina often accompanies disease states associated with suboptimal VMC, including BV, preterm birth, and sexually transmitted infections [51][52][53] . For example, L. crispatus dominance of the vaginal niche is associated with suppressed levels of IL-1β compared to high-diversity compositions 17,54-58 . Further, MBL is a recognition molecule for G. vaginalis 59 and several single-nucleotide polymorphisms in the MBL2 gene have been reported to increase the risk of BV 60,61 . Our findings show that local vaginal immune responses are reflected in the cervicovaginal metabolic phenotype and can be readily detected by DESI-MS. This has clear translational benefit. We have recently reported that cervical cerclage, a procedure used to reinforce the cervix in women at risk of preterm delivery due to cervical shortening, can induce vaginal bacterial dysbiosis, inflammatory activation and premature cervical ripening associated with increased risk of preterm birth, if performed using a braided instead of monofilament cerclage material 22 . Consistent with these findings, DESI-MS showed capacity to detect clinically relevant inflammatory responses to cerclage insertion that associated with patient outcome. Coupled with VMC prediction, this highlights direct swab profiling by DESI-MS as an innovative platform for host-microbiota interactions during pregnancy that could potentially facilitate PTB risk stratification, optimization of preventative interventions 29 (e.g. targeted antibiotic treatment following preterm premature rupture of the foetal membranes and response to treatment interventions during pregnancy 9 ). While our study has focused upon application of DESI-MS in pregnancy cohorts, the rapid detection of VMC and host response offered by this methodology could also inform treatment strategies in other clinical scenarios. For example, efficacy of the topical pre-exposure prophylactic, Tenofovir, is superior in preventing HIV acquisition in women with Lactobacillus-dominated VMCs compared to those dominated by G. vaginalis and BV-associated bacteria 3 . It is therefore conceivable that efficacy and effectiveness could improve through stratification of treatment using point-ofcare metabolic profiling of vaginal swabs. By extension, our approach may also be useful for bedside monitoring of response to treatments designed to optimize VMC, such as live biotherapeutics 62 and vaginal microbiome transplantation 63 .  (IL-1β, IL-8, MBL, C5, C3b/C3bi, IgE, IgG2, IgG3, IgG4, IgM). t-ratio ranges from +10 (red) to −7.5 (blue). C Association between predicted logtransformed value of immune marker by DESI-MS and measured log-transformed values by multiplexed immune-assay for IL-1β (CV R 2 = 0.51), IL-8 (CV R 2 = 0.37), C3b/iC3b (CV R 2 = 0.33), IgG3 (R 2 = 0.31), IgG2 (CV R 2 = 0.27), MBL (CV R 2 = 0.26). A linear regression line was fitted to the logtransformed values and their corresponding prediction. A box plot of predicted immune marker levels for LDEP (red) and LDOM (green) samples is also presented (n = 136 pregnancies, 369 swabs). The lower, interior horizontal line, and upper bounds of the box represent the 25th, median and 75th percentile values, respectively. Whiskers are drawn from the corresponding box boundary to the most extreme data point located within the box bound ± 1.5 × IQR (interquartile range). P values are reported for a two-tailed Welch t-test for the difference in mean predicted immune markers between LDEP and LDOM.
A limitation of our approach was the relatively low sensitivity of prediction. Further inspection of the 16S rRNA gene amplicon data indicated that this was largely associated with misclassification of samples harbouring mixed compositional structures, often containing G. vaginalis (Supplementary Fig. 6). This suggests that hard-clustering techniques often used for determining CSTs may under-estimate the impact of low abundance taxa on the host mucosal metabolome. Comparison of predictive performances from the LC-MS assays indicated that metabolic coverage may also impact misclassification rate with small polar  ). B LDEPL vaginal composition was associated with increased IL-1β levels compared to LDOM; however, highest levels were observed in LDEPL women subsequently having preterm delivery. This relationship was also observed when IL-1β levels and vaginal microbiota composition were predicted using direct swab profiling by DESI-MS (n = 103 pregnancies, 103 swabs). C A relationship between LDEPL, increased MBL and subsequent preterm birth was also detected by DESI-MS profiling (n = 103 pregnancies, 103 swabs). D Elevated MBL and E elevated IL-1β levels were observed in response to cervical cerclage performed with braided cerclage material, but not monofilament material (n = 34 pregnancies, 68 swabs). F Preterm birth in women treated with cervical cerclage using braided cerclage material was associated with higher IL-1β levels compared to term birth outcomes (n = 13 pregnancies, 13 swabs), whereas no relationship between IL-1β levels measured or DESI-MS-predicted were observed with pregnancy outcome following cervical cerclage using monofilament material (n = 21 pregnancies, 21 swabs). All box and whisker plots are drawn with the lower, horizontal interior line, and upper bounds of the box representing the 25th percentile, median and 75th percentile values, and whiskers extending from the lower or upper box bonds to the position of the most extreme data point within ± 1.5 × IQR (interquartile range).  Figs. 5 and 7). It is also worth noting that the analysis of VMC in this study was limited to the assessment of relative abundances of the bacterial component of microbiome and a comparatively small number of patients carrying less prevalent community compositions (e.g., L. gasseri, L. jensenii and B. breve-dominated). The detection of B. brevedominated women in our study was limited to the VMET2 patient cohort where a mixed formulation of the 27F forward primer set was used. This primer formulation has been shown to maintain the rRNA gene ratio of key vaginal species including Lactobacillus spp. to Gardnerella spp. as well as improve detection of Bifidobacterial species, which are otherwise not detected by the 'universal' 27f primer often used in metataxonomics studies 30,64 . Consistent with these results, a number of recent studies have also identified B. breve-dominated vaginal microbiota using metataxonomics approaches 65,66 . Additional datasets that provide improved representation of more diverse compositions and rarer taxa may also eventually facilitate prediction of VMC by DESI-MS beyond major CSTs as described here, and toward species level classification. It is also important to note that recent metataxonomics studies have indicated higher false discovery rates when using relative abundance from NGS data compared to quantitative PCR 67 . Future studies incorporating microbial load within a sample may further improve the sensitivity of DESI-MS predictive models, particularly those with intermediate compositions.
In conclusion, we show that direct on-swab DESI-MS permits rapid and robust assessment of vaginal microbiota:host immune interactions reflected within the cervicovaginal mucosal metabolome. While we show that this may offer an approach for preterm birth risk stratification and selective targeting of preventative treatments, we expect it to have wider application to the assessment of vaginal microbiota-host interactions in other clinical scenarios, including both pregnant and non-pregnant women. Metabolic profiling of cervicovaginal swabs using direct swab analysis by DESI-MS. All chemicals used were analytical reagent grade. HPLC grade methanol and water for DESI-MS analysis and swab sample extraction were purchased from Sigma-Aldrich (St Louis, MO). The method used for the direct analysis of vaginal swab samples by DESI MS is detailed elsewhere 26 . Briefly, we used an LTQ-Orbitrap Discovery mass spectrometer (Thermo Scientific, Bremen, Germany) coupled with a DESI-MS source designed for direct swab analysis. Swabs were placed into a rotating holder positioned orthogonally in front of the MS inlet capillary with a swab-capillary distance of approximately 2 mm. The DESI sprayer tip was pointed to the swab centre with a tip-sample distance of 1.5-2 mm and a distance between the tip and the inlet capillary of 2 mm. The entire surface of the medical swabs was analysed by DESI-MS through clockwise rotation of the swab toward the MS capillary. The cervicovaginal mucosa was absorbed from the swab tip by gently desorbing the analytes with charged droplets of methanol/water (95:5, v/v) mixture and directed to the mass spectrometer. For each sample, 30 scan mass spectra (m/z 50-1000, R = 30,000 (FWHM)) were recorded in the negative and positive ion mode.

Methods
Metabolic profiling of cervicovaginal swabs using direct swab analysis by LC-MS analysis. Liquid extraction was performed on each swab by adding a MeOH:H 2 O (1:1, v-v) solution as eluent to a final concentration of 50 mg vaginal fluid/ml. Each blank swab was extracted with 1 ml solution using a repeated sonication and vortexing cycle for 30 s each. Recovery of soluble material was achieved through centrifugation of swabs (2000 × g for 2 min) seated in a 200 µl loading tip positioned in a sterile microcentrifuge tube. Associated supernatants were pooled and centrifuged at 16,000 × g for 10 min to remove insoluble material before the resulting supernatant was divided into three aliquots. An additional extraction of lipids from the swab was performed by repeating the procedure with an isopropanol:water (4:1, v-v) solution to a final extraction of 25 mg vaginal fluid/ ml. The resulting extract was pooled with one aliquot of the methanol:water extract. All extracts were evaporated using a SpeedVac for further reconstitution before analysis by LC-MS 68,69 .
Reversed phase LC-MS analysis of small metabolites. Sample reconstitution was performed in 300 μl of water. In all, 250 μl aliquot of reconstituted sample material were used for individual sample preparation of 96-well plates including the addition of 25 μl of full RP 2× labelled standard mix (L-glutamine-13 C 5 ; L-glutamic acid 13 C 5 ; creatinine-methyl-D 3 ; cytidine-5,6-D 2 ; citric acid 13 C 6 ; L-isoleucine-13 C 6 15 N; L-leucine-13C 6 ; L-phenylalanine-13 C 9 15 N; hippuric acid-D 5 , benzoic acid-13 C 6 , octanoic acid-13 C 8 , L-tryptophane-13 C 11  ). In addition, 50 μl aliquots of each sample were used for pooling and generation of QC sample. For chromatographic separation a 2 μl aliquot of extracted metabolites from each sample was injected onto a Waters Acquity UPLC BEH C8, 1.7 µm, 2.1 × 100 mm column (Waters Corp.) kept at 55°C using an ACQUITY UPLC system (Waters Corp.). The mobile phase consisting of water:isopropanol:acetonitrile (50:25:25, all high-grade LC-MS solvents from Fisher Scientific or Sigma-Aldrich) + 5 mM ammonium acetate + 0.05% acetic acid + 20 µM phosphoric acid (A) and isopropanol:acetonitrile 50:50 (Sigma-Aldrich) + 5 mM ammonium acetate + 0.05% acetic acid (B). Each sample was resolved for 13.15 min at a flow rate of 0.5 ml/min. Starting conditions were 99% A and 1% B and the gradient changed with a ramp of curve 6 as follows: decrease to 70% A and 30% B over the first 2 min, decrease to 10% A with 90% B from 2 to 11.50 min, decrease to 0.1% A with 99.9% B from 11.50 to 12.50 min, after which the solvent composition returned to starting conditions over 0.1 min until 13.15 min.
HILIC LC-MS analysis. Sample reconstitution was performed in 200 μl of water:acetonitrile mixture (1:3.6, v-v). 115 μl aliquot of reconstituted sample material were used for individual sample preparation of 96-well plates including the addition of 2.5 μl of 48× labelled IS standard mix ((uracil-2-13 C 15 N 2 ; N-benzoyl-d 5glycine, adenosine-2-D 1 , adenine-2-D 1 , taurine-15 N, L-tryptophan-D 5 , L-phenylalanine-13 C 9 15 N, creatine-(methyl-D 3 )-monohydrate, L-arginine-13 C 6 -hydrochloride). In addition, 20 μl aliquots of each sample were used for pooling and generation of QC sample. HILIC chromatography analysis was performed using a Waters Acquity UPLC BEH HILIC (1.7 μm, 2.1 × 150 mm) column (Waters Corporation, Milford, MA, USA) kept at 40°C. The mobile phases consisted of 0.1% formic acid and 20 mM ammonium formate in water (A, Fisher Scientific), and 0.1% formic acid in acetonitrile (B, Sigma-Aldrich) at a flow rate of 0.5 ml/min. Starting conditions were 5% A and 95% B and the gradient changed with a ramp of curve 6 as follows: increase to 20% A and 80% B over the first 5.5 min; increase to 50% A with 50% B from 5.5 to 7.1 min, after which the solvent composition returned to starting conditions over 0.1 min until 12.65 min. The injection volume was 2 μl of extracted metabolites for negative ion mode analysis.
LC-MS instrumental operation and analysis. The column eluent was introduced directly into the mass spectrometer by electrospray. MS was performed on a Waters Xevo G2-S QTOF mass spectrometer (Waters Ltd, Elstree, UK) operating in either negative-ion (ESI−) or positive-ion (ESI+) electrospray ionization mode with a capillary voltage of 1 kV for RP and 1.5 kV for HILIC, and a sampling cone voltage of 20 V. The desolvation gas flow was set to 1000 l/h and the desolvation temperature was set to 600°C. The cone gas flow was 150 l/h, and the source temperature was 120°C. Accurate mass was maintained by introduction of LockSpray interface of Leucine Enkephalin (m/z 556.2771 in ESI+, m/z 554.2615 in ESI−) at a concentration of 200 pg/μl in 50% aqueous acetonitrile with a scan time of 0.07 s over 4 scans, after each interval of 60 s. Data were collected in centroid mode from 50 to 1200m/z in MS scanning mode. To ensure system suitability and stability, a study reference (SR) QC sample was prepared by combining equal aliquots of all the samples and injected at regular intervals throughout the analytical run. This SR sample was also used to condition the column (30 injections) prior to the analysis of both the ESI+ and ESI mode batches. Blank samples (i.e. injection of the reconstitution solvent) were also run to check the presence of artefact or contaminant peaks. A dilution series of the SR sample was also acquired at the beginning and end of injection sequence. Data-dependent acquisition (DDA) and MS E analysis of the SR sample were performed for structural elucidation.
Mass spectral data processing. DESI-MS raw mass spectral data were converted from.raw files into.mzXML format via the ProteoWizard msConverterGUI 70 . The first 30 spectra per sample were averaged into a single spectrum and saved as.csv file using the MALDIquant R package (v1.19.3). The detailed parameter settings and algorithm used in the R package MALDIquant 71 are as follows: Peak detection was performed using 'detectPeaks', with SNR set to 3 and half window size to 10, and using the median absolute deviation (MAD) method. Peaks were aligned with a half window size set to 20, tolerance set to 0.002, signal-to-noise-ratio (SNR) set to 2 using the warping method LOWESS. The reference peaks were created by calling the function 'referencePeaks' from MALDIquant with method 'strict', 'minFrequency' set to 0.01 and peak binning tolerance 0.002.
Metabolite identification. Target m/z features were putatively annotated using online databases including HMDB, Metlin, MMCD, KEGG and Lipidmaps with a 5 ppm tolerance for each compound. To account for potential imprecision in m/z measurement due to peak binning and data processing artefacts, raw m/z value were confirmed by inspection of the.raw data. Metabolite annotation was performed by searching the measured m/z ratios against METLIN (http:// metlin.scripps.edu), Lipidmaps (http://www.lipidmaps.org) and the HMDB (http:// www.hmdb.ca) [74][75][76] . Further structural elucidation was performed using MS/MS experiments via collision-induced dissociation on the LTQ-Discovery MS instrument (Thermo Scientific), and with DDA of the precursor ion by the quadrupole of the Xevo G2 XS Q-TOF mass spectrometer (Waters corporation). For the annotation of metabolites, the MS/MS spectra were matched against spectral libraries from HMDB, NIST and METLIN that were compiled with either authentic standards or theoretical assignment. Identification of metabolites where MS/MS reference spectra were not available were annotated using chemical fragmentation rules. For the structural assignment of glycerophospholipid, fragments of the polar head group or the fatty acyl chains were investigated to confirm the annotation proposed by the databases and discriminate isomers. MSI levels of compound identification were further reported as suggested by the Metabolomics Standards Initiative (MSI) 77 .
Counts from all OTUs assigned to the same species were summed to generate matrices of total counts per species. Species with less than 50 total counts were excluded. Community state types were then assigned to each sample with HC, using Ward-linkage and the Jensen-Shannon distance. The number of clusters which maximized the mean silhouette score were selected. The relative abundance of each species in a cluster was inspected and the clusters matched to CST reported in previous publications, by comparing the relative abundances of the top five OTUs obtained for each cluster. Samples where the most abundant species was a Lactobacillus spp. other than L. crispatus, L. iners, L. gasseri or L. jensenii, where manually assigned to a separate cluster (CST VII). HC analyses and heatmaps were performed using Python (v3.8.5), 'scipy' 83 (v1.5.2), the 'matptlotlib' 84 (v3.3.2) and 'seaborn' 85 (v0.11.0) libraries.
Bacterial isolates were sampled directly from solid agar plates using a medical rayon swab before being transferred into a sterile microcentrifuge tube and stored at −80°C. Three 'blank' swab samples of each solid agar media were also collected. Direct swab analysis using DESI-MS was performed in both positive and negative ion mode. Spectra were processed and assembled into a single data matrix. For each isolate, the intensities per m/z feature of biological replicates were compared against their corresponding culture media background samples using the Wilcoxon-Mann-Whitney test. The list of significant features was matched to the set of identified features from the analysis of human cervicovaginal swab samples. The mean fold change for matched peaks (tolerance < 5 ppm) was estimated as the ratio of the mean intensity within bacterial biomass samples to the mean intensity in the background culture media samples. The heatmap in Supplementary Fig 6 was generated in Python (v3.8.5), with the 'matptlotlib' (v3.3.2) and 'seaborn' (v0.11.0) libraries.
Statistical analysis of metabolomic profiling data. The linear mixed-effects modelling analyses were performed in R using the 'lme4' 87 package (v1. 1.25). For all metabolic profile variables in each data matrix, a linear mixed-effects model with the following 'lme4' formula was fitted: Metabolite~GestationalAge + CST + MaternalAge + BMI + Ethnicity + (GestationalAge||SubjectID). This model structure accounts for the repeated measures by using a random intercept per pregnancy and a random slope per pregnancy. No random effect term was used to model correlation between the random slopes and random intercepts. Gestational age, CST, age, BMI and ethnicity were modelled as fixed effects. All models were fitted using restricted maximum likelihood ('REML = TRUE' in lmer function call). The 'emmeans' 88 package (v1.5.2.1) was used to perform contrast coding and obtain effect size estimates and P values for contrasts and trends of interest. For the comparisons between different CST's and LDOM vs LDEPL, the mean levels and contrasts obtained from the model were estimated assuming maternal age = 30 years, BMI = 23 and gestational age = 20 weeks, and averaged over all levels of ethnicity. For the gestational age trend estimates were further averaged across all CSTs. Detailed tables with effect size estimates and P values for these analyses and contrasts are available via a GitHub repository (see 'Code availability'). The comparisons between Lactobacillus dominant vs depleted was encoded as a comparison between the grand mean of CST I, II III, V and VII levels vs the mean of CST IV. The P values for each contrast tested were calculated using the Kenward-Roger approximation, implemented in the 'pbkrtest' package (v0.4.8.6) 89 . The P values for all contrasts within a single metabolic signal were not corrected for multiple testing. Instead, for each of the main contrasts of interest (early gestation vs late gestation, CST I vs III, CST I vs IV, CST III vs IV and LDOM vs LDEPL), all P values obtained across all metabolic variables in an assay were pooled and FDR corrected together as a signature, using Benjamini-Hochberg false discovery rate correction and selecting a 5% FDR cut-off. The heatmaps in Fig. 1B were created in Python (v3.8.5), 'matptlotlib' (v3.3.2) and 'seaborn' (v0.11.0). Only features which were statistically significant after false discovery rate correction and replicated in both VMET and VMET2 datasets were plotted. To account for differences in m/z calibration between datasets, a feature was only considered to be replicated between datasets if it was possible to find a marker with an m/z error of less than 5 part-per-million in the final FDR corrected signature from the other study. The linear mixed effect model, semi-partial r 2 and conditional R 2 measures were calculated using the R packages 'r2glmm' (v0.1.2) and 'MuMIn' (v1.43.17), respectively. PERMANOVA analyses 90 were performed with the 'vegan' package (v2.5.7), using the adonis2 function, with method = 'euclidean', permutations = 999, and by adding terms sequentially following their order in the model formula log(MetaboliteMatrix + 1)~CST + GestationalAge + Ethnicity + BMI + Age + SubjectID.
For the immune-metabolite association analysis, a series of linear models with the formula lm(log(ImmuneMarker + 1)~Metabolite) were fitted in R (v4.0.3) for all immune-marker cytokine pair combinations, and the corresponding F-test P values calculated. Benjamini-Hochberg false discovery rate correction was used independently for each immune marker analysis, with a 5% FDR cut-off. The heatmaps in Fig. 2A, B were generated using Python (v3.8.5), 'matptlotlib' (v3. 3.2) and 'seaborn' (v0.11.0), using the results of the linear model analysis and the random forest regression models (detailed in the next section).
Prediction of CST, preterm birth and immune markers from the metabolomic data. The CST type was predicted from the metabolomic profiles using random forest classifiers. For the Lactobacillus dominant vs depleted comparison, samples which were assigned to CST I, II, III, V or VII were assigned to the Lactobacillus dominant class, and samples with CST IV or VI label were assigned to Lactobacillus depleted. Random forest classifiers were trained to predict the CST or Lactobacillus depletion status using the metabolic profile variables, without using any other clinical or demographical variables, including gestational age. An assessment of the impact on model performance of using gestational age as covariate was also performed (see the analyses under 'CST Typing by DESI-MS' in the GitHub code repository).
Prediction of preterm birth was also performed using random forest classifiers. The centred-log-ratio transformation was applied to the 16S data matrices with the R package 'propr' (v4.2.6) 91 . The gestational period was divided in three gestational windows (1st timepoint: 0-14 weeks; 2nd timepoint: 14-24 weeks; 3rd timepoint 24-40 weeks) and a random forest classifier models fitted at each window using only the 16S data, immune markers, or DESI-MS profiles as predictors. Only one sample per patient (selected as the closest to the gestational age window mid-point) was used in each model. For the immune marker prediction, immune marker concentrations were initially log-transformed after addition of a constant offset log(x + 1). A random forest regressor was trained to predict a single logtransformed immune marker level. All random forest models were trained in R (v4.0.3), using the 'randomForest' (v4. [6][7][8][9][10][11][12][13][14] 92 package in combination with the 'caret' package (v6.0-86) 93 , for model fitting, performance metric calculation and cross-validation. RF classifiers were trained using constant default parameters: number of trees 'ntree' = 1000, and the number of variables 'mtry = number of metabolic variables/3'. For the random forest regression models, 'mtry = sqrt(number of metabolic variables)'. A repeated (n = 15) stratified fivefold cross-validation procedure was used for all models. This approach ensures that the train and test sets contain a similar proportion of samples from each class.
Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
The metabolic profiling data generated in this study have been deposited and made publicly available in the MetaboLights database under study identifier MTBLS717. The sequence data for the study are publicly available through the European Nucleotide Archive [https://www.ebi.ac.uk/ena] under accession numbers PRJEB 11895, 12577 and (https://www.ebi.ac.uk/ena/browser/view/PRJEB41427). Relevant clinical and patient metadata are publicly available in the GitHub repository at https://www.github.com/ gscorreia89/desims-cst-analysis/, digital resource identifier: https://doi.org/10.5281/ zenodo.5513501 98 .